Wave turbulence and vortices in Bose-Einstein condensation 



Sergey Nazarenko* and Miguel Onorato^ 

* Mathematics Institute, The University of Warwick, 
Coventry, CV4-7AL, UK 
^Dipartimento di Fisica Generale, Universita di Torino, 
Via P. Giuria, 1 - 10125 - 
Torino, Italy 

February 8, 2008 

Abstract 

We report a numerical study of turbulence and Bose-Einstein condensation within the two-dimmensional Gross- 
Pitaevski model with repulsive interaction. In presence of weak forcing localized around some wave number in the 
Fourier space, we observe three qualitatively different evolution stages. At the initial stage a thermodynamic energy 
equipartition spectrum forms at both smaller and larger scales with respect to the forcing scale. This agrees with 
predictions of the the four-wave kinetic equation of the Wave Turbulence (WT) theory. At the second stage, WT 
breaks down at large scales and the interactions become strongly nonlinear. Here, we observe formation of a gas of 
quantum vortices whose number decreases due to an annihilation process helped by the acoustic component. This 
process leads to formation of a coherent-phase Bose-Einstein condensate. After such a coherent-phase condensate 
forms, evolution enters a third stage characterised by three-wave interactions of acoustic waves that can be described 
again using the WT theory. 

1 Background and motivation 

For dilute gases with large energy occupation numbers the Bose-Einstein condensation (BEC) can be described by the 
Gross-Pitaevsky (GP) equation 012: 

j*t + A*-|1'|^v[' = 7, (1) 

where 'I' is the condensate "wave function" (i.e. the c-number part of the boson annihilation field) and 7 is an operator 
which models possible forcing and dissipation mechanisms which will be discussed later. Renewed interest to the 
nonlinear dynamics described by GP equation is related to relatively recent experimental discoveries of BEC 
GP equation also describes light behaviour in media with Kerr nonlinearities. In the nonlinear optics context it is 
usually called the Nonlinear Schoedinger (NLS) equation. 

It is presently understood, in both the nonlinear optics and BEC contexts, that the nonlinear dynamics described 
by GP equation is typically chaotic and often non-equilibrium |51|71|S]|S]. Thus, it is best characterised as "turbulence" 
emphasizing its resemblance to the classical Navier-Stokes (NS) turbulence. On the other hand, the GP model has an 
advantage over NS because it has a weakly nonlinear limit in which the stochastic field evolution can be represented as 
a large set of weakly interactive dispersive waves. A systematic statistical closure is possible for such systems and the 
corresponding theory is called Wave Turbulence (WT) ilOi . For small perturbations about the zero state in the GP 
model, WT closure predicts that the main nonlinear process will be four-wave resonant interaction. This closure was 
used in to describe the initial stage of BEC. In the present paper we will examine this description numerically. 

We report that our numerics agree with the predicted by WT spectra at the initial evolution stage. 



It was also theoretically predicted that the four-wave WT closure will eventually fail due to emergence of a coherent 
condensate state which is uniform in space At a later stage the condensate is so strong that the nonlinear dynamics 
can be represented as interactions of small perturbations about the condensate state. Once again, one can use WT 
to describe such a system, but now the leading process will be a three-wave interaction of acoustic-like waves on the 
condensate background W. Coupling of such acoustic turbulence to the condensate was considered in 111) which allowed 
to derived the asymptotic law of the condensate growth. However, this picture relies on assumptions that the system 
will consist of a uniform condensate and small perturbations. Neither the condensate uniformity nor the smallness of 
perturbations have ever been validated before. In the present paper we will examine whether it is true that the late 
stage of GP evolution can be represented as a system of weakly nonlinear acoustic waves about strong quasi-uniform 
condensate. By examining the frequency-wave number Fourier transforms, we do observe waves with frequency in 
agreement with the Bogoliubov dispersion relation. The width of the frequency spectrum is narrow enough for these 
waves to be called weakly nonlinear. 

A big unresolved question in the theory of GP turbulence has remained about the stage of transition from the four- 
wave to the tree-wave regimes. This stage is strongly nonlinear and, therefore, cannot be described by WT. However, 
using direct numerical simulations of equation Q, we show that the transitional state involves a gas of annihilating 
vortices. When the number of vortices reduces so that the mean distance between the vortices becomes greater than 
the vortex core radius (healing length) the dynamics becomes strongly nonlinear. This corresponds to entering the 
Thomas-Fermi regime when the mean nonlinearity is greater than the dispersive term in the GP equation. The mean 
inter-vortex distance is a measure of the correlation length of the phase of and, therefore, the vortex annihilation 
corresponds to creation of a coherent-phase condensate. At this point, excitations with wavelengths in between of the 
vortex-core radius and the inter-vortex distance behave as sound. In this paper, we draw attention to the similarity of 
this transition process to the Kibble-Zurek mechanism of the second-order phase transition in cosmology |12II13| . 

2 WT closure and predictions 

The WT closure is based on the assumptions of small nonlinearity and of random phase and amplitude variables. To 
show how these assumptions come about we reproduce several essential steps of the kinetic equation derivation. Let 
us consider GP equation Q in a periodic box and write it in the Fourier space: 

idt'^k - = J2 *"*M*^<5;;^ + 7fc, (2) 

where "i/j = 5'(kj), overbar means complex conjugation, wave vectors kj {j — 1,2,3) are on a 2D grid (due to 
periodicity) and symbol S^" = 1 for k -I- = k^ -|- ki^ and equal to otherwise. 

2.1 Four- wave interaction regime 

In order to describe the WT theory for equation (|5J we will neglect the forcing/dissipation term 7^ assuming that 
those are localized at high or low wave numbers and we are mainly interested in an inertial range of k. Let us make a 
transformation to the interaction representation variables at which absorb the linear dynamics and oscillations due to 
the diagonal terms k = k2 and ki = ka, 

flfe =6""''*+"''""'**^, (3) 

where LUk ~ is the linear frequency and 

k 

is the nonlinear frequency shift. In terms of ak we have 

idtUk ^ ^2 raaaMai- (5^" e'""" (5) 

where co^" = ujk + (j->a ^ <-Ofj, ~ <-l>v ■ The nonlinear frequency drops out from this expression because it is fc- independent. 
Small nonlinearity allows to separate the linear and nonlinear timescales. Let us consider solution a^iT) at an inter- 
mediate (between linear and nonlinear) time T ^ 2Ti/ujk but such that changes in at are still small. We seek this 
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solution as a series 

O-k — O,^, + Clf. + tlj, + ••■ (.oj 
which is obtained by recursive substitution into 0. For example, a^^^ — a^lr^o and 

at\T) = ^^Y^-a^^^a^^\'^:^^5lZAlZ. (7) 
where Ajj^ = (e'"^*"'"^ — Vj/iuo^"^. For the second iteration we get 

(2)/Tn\ / xM" ^feu (0) (0) (0)-(0)-(0) ,^1/ k^lv feux nrau rfeii-(0)-(0) (0) (0) (0)^1/ kau ku\\ /on 

where E{x,y) = J^^ A(a; — y)e^^^dt. Expressions Q and Q are sufficient for writing the leading order (in small 
nonlinearity) for the evolution of the spectrum, — (jafej^): 

Uk = {{\ak{T)\') - (|afe(0)|'))/T = ^i{af^a[''>) + (ai'^af) + {af^ai'^} + (af af) + {|aW|^)). (9) 

Here, we must substitute expressions Q and Q and perform averaging over the ensemble of initial fields a]^^ ■ For this, 
we will use a generalized RPA (Random Phase and Amplitude) approach introduced in |i4) and which is different from 
the traditional RPA (Random Phase approximation, see e.g. |10p by emphasizing randomness of the amplitudes and 
not only the phases. Namely, we represent complex amplitudes as a^"' = Akipk where Ak are positive real amplitudes 
and V'fe are phase factors taking values on the unit circle on the complex plane. Generalized RPA says that the initial 
wave field is such that all of its amplitudes Ak and phase factors ipk make a set of independent random variables and 
that ipk's are uniformly distributed on the unit complex circle (distribution of Ak's need not be specified). 

RPA allows to close equations for the spectrum JSJ by using the Wick-type splitting of the higher Fourier moments 
in terms of the spectrum. The first two terms in @ contain three a*^"' 's each and, therefore, vanish due to the phase 
randomness. The other three terms, after substitution of Q and 101, RPA averaging and taking the large-box and 
large T limits, give (see details e.g. in |14p: 

nfc = 47r / nkUun^n,, (— - + — " — ) (5(<jJa")<5M" dkudk^dk^. (10) 

J \nk Uu J 

This is the wave-kinetic equation (WKE) which is the most important object in the wave turbulence theory (for the 
GE equation, it was first derived in 0). It contains Delta functions for four wave vectors, 5^" — 5(k + k^ — k^j — k,y), 
and for the four corresponding frequencies, 5(tjp"), which means that the spectrum evolution in this case is driven 
by a 4-wave resonance process. Note that WT approach is applicable not only to the spectra but also to the higher 
moments and even the probability density functions |14l I15| . However, we are not going to reproduce these results 
because their study is beyond the aims of the present paper. 

There are four power-law solutions of the 4-wave kinetic equation in llUl and they are related to the two invariants 
for such systems, the total energy, E = J ujkrikdk, and the total number of particles, N — f n^dk. Two of such 
power-law solutions correspond to a thermodynamic equipartition of one of these invariants, 

Uk ~ l/'i'fe = (energy equipartition), (11) 

rik = const (particle equipartition). (12) 

These two solutions are limiting cases of the general thermodynamic distribution, 

Uk = T/{ujk + (13) 

where constants T and /i have meanings of temperature and chemical potential respectively. Due to isotropy, it is 
convenient to deal with an angle-averaged ID wave action density in variable k = |k|, so called ID wave action 
spectrum A'^fe = 2nknk. In terms of Nk, solutions IIH and 1121 have exponents —1 and 1 respectively. 

The other two power-law solutions correspond to a Kolmogorov-like constant fiux of either energy (down-scale 
cascade) or the particles (up-scale cascade) As shown in 0, the formal solution for the inverse cascade has wrong 
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sign of the particle flux and is, therefore, irrelevant. On the other hand, the power exponent of the direct cascade 
solution formally coincides with the energy equipartition exponent —2 and, in fact, it is the same solution. Because of 
such a coincidence, the energy flux value is equal to zero on such a solution and, therefore, it is more appropriate to 
associate it with thermodynamic equilibrium rather than a cascade. 



2.2 Three-wave interaction regime 

If the system is forced at large wave numbers and there is no dissipation at low k's then there will be condensation 
of particles at large scales. The condensate growth will eventually lead to a breakdown of the weak nonlinearity 
assumption |5] 117) and the 4-wave WKE UlUl will become invalid for describing subsequent evolution. On the other 
hand, it was argued in that such late evolution one can consider small disturbances of coherent condensate state 
\l/o ~ const, so that a WT approach can be used again (but now on a finite-amphtude background), 

*(x,t) = *o(l + <?i(x,t)), ?i«*o. (14) 

Then, with respect to condensate perturbations 0, the linear dynamics has to be diagonalised via the Bogoljubov 
transformation, which in our case is |^ 1161 ITT] 



1/2 \ / ,.,1/2 

k 



1/2 ' fc ' ,,,1/2 k 



(15) 



where are new normal amplitudes (see for example lllp. In the linear approximation, amplitudes Ofe oscillate at 

frequency 

ujk = k^WTWo (16) 
which is called the Bogoljubov dispersion relation. For strong condensate, po ^ fc^, this dispersion relation corresponds 
to sound. 

Because of the non-zero background, the nonlinearity will be quadratic with respect to the condensate perturbations 
and, thus, the resulting WT closure now gives rise to a 3-wave WKE. This WKE was first obtained in [2| (see also 
[lip and here we reproduce it without derivation. 



where 



TT j {Rki2 - Rik2 - R2ki) dkidka , (17) 
-Rfci2 = |Vfcfcifc2 12 (5(k — ki — k2) S{ujk —ioi— UJ2) {n\n2 ~ n^ni — 71^712). 



VkMM - (27r)3/2 \(a,Qia2)i/2 2 



(k-kQ ^ (k-k,) ^ (ki-k,) 



kk\oi2 fcfeai k\k2Cik 



(18) 



where 



aft = 2po + fc^ (19) 

At late time the condensate becomes strong, pQ 3> fc^, and turbulence becomes of acoustic type. The number of 
particles is not conserved by the turbulence alone (particles can be transferred to the condensate) and there are only 
two relevant power-law solutions in this case: thermodynamic equipartition of energy and the energy cascade spectrum. 
Because of isotropy, one often considers ID (i.e. angle-integrated) energy density, 

E{k) = 2nkuJknk. (20) 

In terms of this quantity, the thermodynamic spectrum is 

E{k) ~ k, (21) 

and the energy cascade spectrum is 

E{k) ~ fc-^/^ (22) 



4 



Note that the energy cascade is direct and the corresponding spectrum can be expected in fe's higher than the forcing 
wave number, whereas the thermodynamic spectrum is expected at the low-fc range to the left of the forcing 
Note that the described above picture of acoustic WT rehes on two major assumptions. 

1. Condensate is coherent enough so that its spatial variations are slow and it can be treated as uniform when 
evolution of the perturbations about the condensate is considered. In the other words, a scale separation between the 
condensate and the perturbations occurs. 

2. Coherent condensate is much stronger than the chaotic acoustic disturbances. This allows to treat nonlinearity 
of the perturbations around the condensate as small. 

Both of these assumptions have not been validated before and their numerical check will be one of our goals. 
Another major goal will be to study the transition stage that lies in between of the 4-wave and the 3-wave turbulence 
regimes. This transition is characterised by strong nonlinearity and the role of numerical simulations becomes crucially 
important in finding its mechanisms. 

Once the 3-wave acoustic regime has been reached, the condensate continues to grow due to a continuing influx of 
particles from the acoustic turbulence to the condensate. This evolution, where an unsteady condensate is coupled with 
acoustic WT, was described in who predicted that asymptotically the condensate grows as po ~ if the forcing is 
of an instability type 7 = Vk^k- However, in the present paper we work with a different kind of forcing which is most 
convenient and widely used in numerical simulations: we keep amplitudes in the forcing range fixed (and we chose their 
phases randomly). Thus, one should not expect observing the regime predicted in in our simulations. Note that 
2D NLS turbulence was simulated numerically with specific focus on the condensate growth rate in In our work, 
we do not aim to study the condensate growth rate because it is strongly dependent on the forcing type which, in our 
model, is quite different from turbulence sources in laboratory. On the other hand, we believe that the main stages of 
the condensation, i.e. transition from a 4-wave process, through vortex annihilations, to 3-wave acoustic turbulence, 
are robust under a wide range of forcing types. 

3 Setup for numerical experiments 

In this paper we consider a setup corresponding to homogeneous turbulence and, therefore, we ignore finite-size effects 
due to magnetic trapping in BEC or to the finite beam radii in optical experiments. For numerical simulations, we 
have used a standard pseudo-spectral method for the 2D equation The number of grid points in physical space was 
set to N X N with N = 256. Resolution in Fourier space was Afc — 2-k/N. Sink at high wave numbers was provided by 
adding to the right hand side of equation Q the hyper- viscosity term iy{—\7'^)"ip- Values of u and n were selected in 
order to localized as much as possible dissipation to high wave numbers but avoiding at the same time the bottleneck 
effect. We have found, after a number of trials, that = 2 X 10"^ and n = 8 were good choices for our purposes. In 
some simulations, we have also used a dissipation at low wave numbers of the form of V^)~"^ with u — 1 x 10~^* 
and n — 8. This was done, e.g., to see what changes if one suppresses the condensate formation. Forcing was localized 
in Fourier space and was chosen as / = |/|ea;p[— with |/| constant in time and </>(t) randomly selected between 
and 27r each time step, i) To study turbulence in the down-scale inertial range we force the system isotropically at wave 
numbers 4Afc < < 6AA:. To avoid condensation at large scales we introduce a dissipation at low wave numbers, as 
was previously explained, ii) To study the condensation we chose forcing at wave numbers 60Afc < \k\ < 63Afc and 
dissipation at all higher wave numbers. A number of numerical simulations were performed both with and without 
the dissipation at the low wave numbers. Time step for integration was t — 0.1 and usually 10^ time steps have been 
performed for each simulation. This is usually enough for reaching a steady state when dissipation at both high and 
low wave numbers was placed. 

4 Numerical results 

4.1 Turbulence with suppressed condensation 

We start with a state without condensate for which WT predicts four-wave interactions. WKE has two conserved 
quantities in this case, the energy and the particles, and the directions of their transfer in the scale space must be 
opposite to each other. Indeed, let us assume that energy flows up-scale and that it gets dissipated at a scale much 
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greater than the forcing scale. This would imply dissipating the number of particles which is much greater than what 
was generated at the source (because of the factor difference between the energy and the particle spectral densities). 
This is impossible in steady state and, therefore, energy has to be dissipated at smaller (than forcing) scales. On the 
other hand, the particles have to be transferred to larger scales because dissipating them at very small scales would 
imply dissipating more energy than produced by forcing. This speculation is standard for the systems with two positive 
quadratic invariants, e.g. 2D Euler turbulence where one invariant, the energy, flows up-scale and another one, the 
enstrophy, flows down-scale. 

Thus, ideally, one would like to place forcing at an intermediate scale and have two inertial ranges, up-scale and 
down-scale of the source. However, this setup is unrealistic because the presently available computing power would 
not allow us to achieve simultaneously two inertial ranges wide enough to study scaling exponents. Therefore, we split 
this problem to two, with forcing at the left and at the right ends of a single inertial range. 

4.1.1 Turbulence down-scale of the forcing 

Our first numerical experiment is designed to test the WT predictions about the turbulent state corresponding to the 
down-scale range with respect to the forcing scale. Thus we chose to force turbulence at large scales and to dissipate 
it at the small scales as described in the previous section. Our results for the one-dimensional wave action spectrum is 
statistically stationary condition is shown in Figure^ We see a range with slope —1 predicted by both the Kolmogorov- 
Zakharov (KZ) energy cascade and the thermodynamic energy equipartition solutions of the four-wave WKE. As we 
mentioned earlier, it would be more appropriate to interpret this spectrum as a quasi-thermodynamic state rather than 
the KZ cascade because the energy flux expression formally turns into zero at the power spectrum with —1 exponent. 
We emphasise, however, that the state here is quasi-thermodynamic with a small flux component present on thermal 
background because of the presence of the source and sink. One could compare this state to a lake with two rivers 
bringing the water in and out of the lake. In comparison, pure KZ cascade would be more similar to a waterfall. 

To check that the waves in our system are indeed weakly nonlinear, we look at the space-time Fourier transform 
of the wave field. The frequency- wave number plot of this Fourier transform is shown in Figure |21 We see that this 
Fourier transform is narrowly concentrated near the linear dispersion curve, which confirms that the wave field is 
weakly nonlinear. We can also see that the spectrum is slightly shifted upwards by a value which agrees with the 
nonlinear frequency shift found via substitution of the numerically obtained spectrum into 

4.1.2 Up-scale turbulence 

In the up-scale range one could expect that, in analogy with the 2D Navier-Stokes turbulence, there would be an inverse 
cascade of the number of particles and that the corresponding KZ spectrum would be observed. Nevertheless, it was 
pointed out in |3| that the analytical KZ spectrum has "wrong" direction of the flux of particles in the 2D GP model 
and, therefore, cannot form. Our numerics agree with this view. Instead of the KZ, our numerical simulations show 
that a statistical stationary state with a power law very close to forms, see Figure |3 This solution corresponds 
to the thermodynamics solution with energy equipartition in the fe-space. Note that both theoretical rejection of the 
particle-cascade spectrum |^ and our numerical study relate to the 2D model and the situation can change in the 3D 
case. ^ Namely, it is possible that the up-scale dynamics in 3D will be characterised by the particle-flux KZ solution 
or a more complicated mixed state which involves both cascade and temperature. On the other hand, formation of a 
pure thermodynamic state in 2D is quite fortunate for the theoretical description because analogies with the theories 
of phase transition between different types of thermodynamic equilibria become more meaningful. 

Here, we also check that the waves in this regime are weakly nonlinear by looking at the space-time Fourier 
transform. The corresponding frequency- wave number plot is shown in Figure 2] As in the down-scale inertial range, 
we see that this Fourier transform is narrowly concentrated near the linear dispersion curve, i.e. the wave field is 
weakly nonlinear in this state. 

^Another difference with the 3D case may be that in 3D the condensate forms as a sharp peak at the ground state whereas in 2D no such 
peak is observed 1191 . Such absence of a sharp peak at fc = in 2D is consistent with our numerics. 
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Figure 3: ID waveaction spectrum Nk in the up-scale range. A power law of the form of A; ^ is also shown. 




Figure 4: Wave number-frequency distribution of the space-time Fourier transform of 4* in the up-scale inertial range. 
Dispersion relation from linear theory is shown in black curve. 



8 




Figure 5: Initial stages of the evolution of the ID wave action spectrum Nk- A power law of the form of A; ^ is also 
shown. 

4.2 Bose-Einstein condensation 

4.2.1 Initial stage: four-wave process 

In order to study the stages of the condensation process, the results presented in the following have been obtained with 
forcing localized at high wave numbers without dissipation at low wave numbers. At the initial stage of the simulation, 
the nonlinearity remains small compared to the dispersion in the GP equation and the four-wave kinetic equation can 
be used. In Figure |3 we show the initial (pre-condensate) stages of the spectrum evolution. Similarly to the case 
where the condensation was suppressed, we observe the formation of a thermodynamic distribution. 

4.2.2 Transition 

After the stage where the four-wave interaction dynamics holds, the dynamics is characterised by a transitional stage 
in which the the low-fc front of the evolving spectrum reaches the largest scale (at about t = 5000), see Figure [H] 
the spectrum begins to become steeper at low wave numbers and, as expected, the thermodynamics solution does 
not hold anymore. This behaviour indicates that a change of regime occurs around time t = 5000. However, the 
information contained in the spectrum is insufficient to fully characterize this regime change and this brings us to 
study this phenomenon by measuring several other important quantities. 

To get an initial impression of what is happening during the transition stage it is worth to first of all examine the 
field distributions in the coordinate space. Figure |7] shows a series of frames of the real part of 'I' (imaginary part looks 
similar). One can see that this field exhibits growth of large-scale structure. 

On the other hand, field |^|, shown in Figure|H| still remains dominated by small-scale structure. In contrast with 
I*!*!, field contains an additional information - the phase. Thus, separation of the characteristic scales in Figure |7| 
andEI can be attributed to the fact that the phase correlation length becomes much longer than the typical wavelength 
of sound (characterised by fluctuations of l^*] as explained above in Section r2.21 . This scale separation can also be seen 
by comparing the spectrum of l^*], shown in Figure |3] with the spectrum of 'i' in Figures |S] and [H] one can see that the 
former is more flat than the later. Now that we have established that the phase is an important parameter, we can 
measure its correlation length as the mean distance between the phase defects - vortices. Vortices in the GP model 
are points in which ^' = 0. Some of such points correspond to the 2n phase increment when one goes once around 
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Figure 6: Later stages of the evolution of the ID wave action spectrum Nk- A power law of the form of fc is also 
shown. 
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Figure 8: \ij{x, y)\ at different times: t = 2500, t = 5000, t = 7500, t = 10000. 




Figure 9: Spectrum for variable |'!/'(a:,y)| at different times. 
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them, whereas the other points gain —2tt. These vortices can be defined as positive and negative correspondingly. In 
contrast with the Euler equation of the classical fluid, positive and negative vortices can annihilate in the GP model 
and they can get created "from nothing" . Figure 1101 shows a sequence of plates showing the positive and negative 
vortex positions at several different moments of time. One can see that initially there were a lot of vortices, which is 
not surprising because initial field is weak, i.e. close to zero everywhere. However, at later times we see the number of 
vortices is rapidly dropping, which means that the vortex annihilation process dominates over the vortex-pair creations. 
The total number of vortices (normalised by A'^^) is shown as a function of time in Figure [TTI where one can see a fast 
decay. The law of decay is best seen on the log — lin plot, see Fiigure [T?l where one can see a regime 

N-uortices = A~ Blogt (23) 

with A = 3.36 and B — 0.9223 which sets in at t = 800 to t = 3500.^ Thus, the phase correlation distance, being of 
the order of the mean distance between the vortices, exhibits a fast growth in time. 

Let us have a look at a slice of the field l^*] through typical vortices at a late time when most of the them have 
annihilated, see Figure [T^ One can see that l^*] is close to zero (i.e. both Re[ilj\ and /m[7/;] cross zero) at the vortex 
centres and that it sharply grows to order-one values ("heals") at small distances from the vortex centres which are 
much less than the distance between the vortices. This means that these vortices represent fully nonlinear coherent 
structures each of which can be approximately seen as an isolated Pitaevski vortex solution j^. In contrast, the initial 
vortices are too close to each other to be coherent and they correspond to a nearly linear field. The moment when 
the mean inter-vortex separation becomes comparable to the healing length can be captured by the intersection point 
of the graphs for the mean (space averaged) nonlinear and the mean (space averaged) Laplacian terms in the GP 
equation, see Figure WW This intersection (at t = 6950) marks the moment when mean nonlinearity becomes greater 
than the mean linear dispersion, i.e. the Thomas-Fermi regime sets in. This regime could be thought of as the one of a 
fully developed condensate when the nonlinearity, when measured with respect to the zero level, is strong and therefore 
the 4-wave WT description breaks down. However, as we will see in the next section, we now have weakly nonlinear 
perturbations if they are measured with respect to a non-zero condensate state. Evolution of such perturbations takes 
form of three-wave acoustic turbulence. 

What makes vortices annihilate? A positive-negative vortex pair, when taken in isolation, would propagate with 
constant speed without changing the distance between the vortices |23| . Thus, there should be an additional entity 
which could exchange energy and momentum with the vortex pair and to allow them annihilate. We note that the 
field I*!'! is very "choppy" in the region between the vortices, see figure 113II . and, therefore, it is natural to conjecture 
that the missing entity is sound. To check this conjecture, we perform the following numerical experiment. At a 
desired time we filter the field and let it evolve further without sound. The filtering is performed numerically in the 
following way: we have used a Gaussian filter in physical space and have smoothed the field around vortices. The 
complex field V is therefore convoluted with a normalised Gaussian function with standard deviation much smaller 
with respect to the mean distance between vortices. The filter is applied only in the region where no vortices are 
located. The result of the filtering procedure on the evolution of the number of vortices is shown in Figure IT^ We 
see that removing the sound component does indeed reverse the vortex annihilation process and for some time (until 
new sound gets generated from forcing) we observe that the vortex creation process dominates. We point out that the 
described above regime change, accompanied by vortex annihilations, is very similar to the Kibble-Zurek mechanism of 
the second-order phase transition in cosmology fl2lll3| . This mechanism suggest that at an early inflation stage, Higgs 
fields experience a symmetry breaking transition from "false" to "true" vacuum, and this transition is accompanied by 
a reconnection-annihilation of "cosmic strings" which are 3D analogs of the 2D point vortices considered in this paper. 
To describe these fields, one normally uses nonlinear equations of so-called Abbelian model |21|. but the non-linear 
Klein-Gordon or even the GP equation are sometimes used as simple models in cosmology which retain similar physics 

4.2.3 Late condensation stage: acoustic turbulence 

It was predicted in ^ that the turbulent condensation in the GP model will lead to creation of a strong coherent mode 
with fc = such that the excitations at higher wave numbers would be weak compared to this mode. If this is the case, 

^At present, we do not have a theoretical explanation of this law of decay. 
^For this reason such vortices are sometimes called "ghost vortices" 1201 . 
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Figure 10: Vortices in the {x, y) plane at different times: t = 2500, t = 3250, t = 5000, t = 7500. 
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Figure 11: Evolution in time of the density of vortices in a lin-log plot. 
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Figure 12: Evolution in time of the density of vortices in a log-lin plot. The dashed hne corresponds to the fit Nyorti, 
3.36 - 9.223Log{t) 
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Figure 13: Slice of the field j^*] for constant y: a single vortex is visible in the plot. 




Figure 14: The solid line represents the space-averaged \V'^tp{x,y)\; the dotted line is the space-averaged \il){x,y)\^. See 
text for comments. 
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Figure 15: Evolution of the vortex density in time. At time t = 6500, sound has been filtered according to the methodology 
described in the text. 



one can expand the GP equation about the new equilibrium state, uniform condensate, use Bogoljubov transform to 
find new normal modes and a dispersion relation for them, equation H15|l . and to obtain a new WKE for this system 
that would be characterised by by three- wave interactions, 1171 1. However, as we saw in Figure |S| the peak at small 
k remains quite broad that is the coherent condensate, if present, remains somewhat non-uniform. Despite of this 
non-uniformity, one can still use the approach of if there is a scale separation between the condensate coherence 
length (intervortex distance) and the sound wavelength and if the sound amplitude is much smaller than the one of 
the condensate. 

We have already seen tendency to the scale separation in figures|S|[7||H||n| On the other hand, smallness of the sound 
intensity can be seen in figure IT^ which compares (space-averaged) < l^*]^ > and < I"!'] >^. We see that at the late 
stages these quantities have very close values which means that the deviations of from its mean value (condensate) 
are weak. Thus, both conditions for the weak acoustic turbulence to exist are satisfied at the late stages. However, the 
best way to check if the condensate perturbations do behave like weakly nonlinear sound waves obeying the Bogoljubov 
dispersion relation consists in plotting the square of the absolute value of the space-time Fourier transform of 'P. This 
result is given in Figure [T71 for the latest stage of the simulation (from time t — 10488 to t = 11000). Note that for 
each k the spectrum has been divided by its maximum in order to be able to follow the dispersion relation up to high 
wave numbers. 

The normal variable for the Bogoljubov sound is given in terms of tp by expressions (1141 and 11511 , and, therefore, 
when plotting the Bogoljubov dispersion 1161 . we should add a constant frequency of the condensate oscillations, 
ujo ~ (l^'l^). One can see that the main branch of the spectrum does follow the Bogoljubov law up to the wave 
numbers which correspond to the dissipation range.** Further, the wave distribution is quite narrowly concentrated 
around the Bogoljubov curve which indicates that these waves are weakly nonlinear. Note that the lower branch in 
Figure IT7I is related to the contribution to expression 11511 which vanishes at larger k. Importanlty, we can also see 
the middle (horizontal) branch with frequency ujo which quickly fades away at finite fc's and which corresponds to the 
coherent large-scale condensate component. 

Now let us consider the energy spectrum. The GP Hamiltonian can be written in terms of both real and Fourier 

*At the same time, these wave numbers are of order of the inverse healing length, and it is unclear whether the Bogoljubov mode is not 
seen there due to the wave dissipation or due to contamination of this range by the broadband (in frequency) vortex motions. 
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Figure 18: Ei{k) at the latest stage of the snTiulation (see equation 12411 '). 

quantities, 

H = l (|V*r + dx = I {k'W' + ^\pf) dk, (24) 

where p = l^'p. Thus, we measure the ID energy spectrum in this case as E{k) = E\(k) + E2{k) with Ei{k) — 

and E2{k) = ||/3p The contributions to the energy spectrum Ei and E2 as well as the total spectrum E{k) at a time 

corresponding to the acoustic regime are shown in figures 

EHl m and 1201 respectively. 

We see that at small scales the total energy spectrum E{k) scales as 1/fc which is a thermodynamic energy equipar- 
tition solution in this case. 

5 Conclusions 

Firstly, we confirmed WT predictions of the energy spectra in the down-scale and up-scale inertial intervals in the cases 
when the fluxes are absorbed by dissipation in the end of the inertial interval (so that no condensation or build-up is 
happening). In both of these cases we observed spectra with an exponent corresponding to the energy equipartition 
thermodynamic solution Nk ~ (which formally coincides with the exponent for the energy cascade solution). By 
looking at the shape of the frequency-wave number mode distributions, we verified that the turbulence is weak. 

Secondly, we studied a system without dissipation at large scales. We observed a process of Bose-Einstein conden- 
sation and formation of a coherent large-scale mode which happens via annihilating vortices. This scenario is similar 
to the Kibble-Zurek phase transition in cosmology which involves annihilating cosmic strings. We established that the 
process of the vortex annihilation is due to the presence of sound. The condensate correlation length, which in our 
case is of the order of the mean inter- vortex distance, turns into infinity in a finite time as A ~ l/(logt* — logt)^''^, see 
equation I23I I. 

We confirmed numerically that in late condensation stages the system can be described as a weakly nonlinear 
acoustic turbulence on the background of a quasi-uniform coherent condensate. Namely, we confirmed that the wave 
excitations are narrowly distributed around the Bogoljubov dispersion law, i.e. that the turbulence is (i) acoustic and 
(ii) weak. We observed a spectrum that corresponds to the energy equipartition solution of the 3-wave kinetic equation 
for such acoustic turbulence. 

An interesting question to be addressed in future is to what extend the findings of this work are relevant to the 
3D GP model. We can speculate that the energy spectra may have a different nature in 3D and, in particular, may 
expect formation of the Kolmogorov-like spectra corresponding to the energy and the particle cascades. On the other 
hand, it is reasonable to expect that the Kibble-Zurek scenario of condensation will persist in the 3D case, i.e. the 
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correlation length will grow because of the reconnecting and shrinking vortex loops. It is also likely that such vortex 
loop shrinking will be facilitated by the sound component. Computations of 3D GP equation in a non-turbulent setting 
were done in |24| where such processes as vortex reconnection and the role of the acoustic component were considered. 
Turbulent setting will be more taxing on the computing resources due to the great variety of scales involved and, 
therefore, necessity of high resolution and long computation times. 
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